TUM-HEP-524/03 



Variations on KamLAND: likelihood analysis 
and frequentist confidence regions 



Thomas Schwetz 



Institut fur Theoretische Physik, Physik Department, Technische Universitat 
Miinchen, James-Franck-Strasse, D-85748 Garching, Germany 



Abstract 



In this letter the robustness of the first results from the KamLAND reactor neutrino 
experiment with respect to variations in the statistical analysis is considered. It is 
shown that an event-by-event based likelihood analysis provides a more powerful 
tool to extract information from the currently available data sample than a least- 
squares method based on energy binned data. Furthermore, a frequentist analysis of 
KamLAND data is performed. Confidence regions with correct coverage in the plane 
of the oscillation parameters are calculated by means of a Monte Carlo simulation. I 
find that the results of the usually adopted x^-cut approximation are in reasonable 
agreement with the exact confidence regions, however, quantitative differences are 
detected. Finally, although the current data is consistent with an energy independent 
flux suppression, a ~ 2a indication in favour of oscillations can be stated, implying 
quantum mechanical interference over distances of the order of 200 km. 
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1 Introduction 



The outstanding results from the KamLAND reactor neutrino experiment 
have lead to a significant progress in neutrino physics. The observed disap- 
pearance of reactor anti-neutrinos is in agreement with the so-called LMA 
solution of the solar neutrino problem |2| . Alternative oscillation solutions like 
LOW, VAC or SMA are ruled out with very high confidence level jS HI IS IS 

IS BE HQ, and 

non-oscillation mechanisms can play only a sub-leading 
role (for a review and references see Ref. 
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These important conclusions are based on a data sample consisting of 54 events 
above the geo-neutrino threshold in KamLAND. The purpose of this letter is 
to discuss issues related to the statistical analysis of these data. In Sec. 2 an 
event-by-event based likelihood analysis is compared to the widely used least- 
squares method based on energy binned data. It is shown that the likelihood 
method allows one to extract more precise information about the oscillation 
parameters from KamLAND data. Since the currently available data sam- 
ple consists only of rather few events, one might ask the question whether the 
approximate confidence regions obtained from the usual ^-cvX method are re- 
liable. In Sec. 3 this question is addressed by calculating frequentist confidence 
regions for the oscillation parameters according to the prescription given by 



Feldman and Cousins [14j . The explicit construction of the confidence regions 
by Monte Carlo simulation takes properly into account statistical fluctuations 
of the rather small data sample and the non-linear character of the oscillation 
parameters. In Sec. 4 the statistical significance of an oscillatory signal in the 
KamLAND data is discussed, and I conclude in Sec. 5. 



2 Comparing likelihood and least-squares methods 



The current KamLAND data sample consists of 86 anti-neutrino events in the 
full energy range. In the lower part of the spectrum there is a relevant contri- 
bution from geo-neutrino events to the signal. To avoid large uncertainties as- 
sociated with the geo-neutrino flux an energy cut at 2.6 MeV prompt energy is 
applied for the oscillation analysis, and 54 anti-neutrino events remain in the fi- 
nal sample. All analyses of KamLAND data SBH00BBI3QHII1II3 



performed so far outside the KamLAND collaboration are using these data 
binned into 13 energy intervals above the geo-neutrino cut, as given in Fig. 5 
of Ref. 0. 1 In Subsec. 2.1 I describe an alternative analysis based on the like- 
lihood function of the data, which allows one to take into account the precise 
energy information contained in each single event. The results of this analy- 
sis are compared to the ones from the energy binned least-squares method in 
Subsec. 2.2. 

Before exact confidence regions are calculated in Sec. 3 the usual x 2_cu t ap- 
proximation will be used. One constructs a statistic AX 2 (sin 2 29, Am 2 ) from 
the data. Under certain assumptions, like large sample limit and linear pa- 
rameter dependence, this statistic will be distributed as a \ 2 with 2 degrees 
of freedom, independent of the point in the parameter space Then a 



For an analysis including the geo-neutrino events see Ref. [1 
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given point (sin 2 29, Am 2 ) is contained in the allowed region at ft CL if 

AX 2 (sm 2 29,Am 2 ) <A X l(2), where / fy(x, n)dx = ft . (1) 

<J 

Here f x ^(x, n) denotes the ^-distribution with n degrees of freedom. In the 
following I will refer to this procedure as "x 2 -cut method". In this section 
it will be applied to calculate approximate confidence regions by using the 
likelihood as well as the least-squares method. 



2.1 Likelihood analysis of KamLAND data 



For given oscillation parameters sin 2 29 and Am 2 the predicted event spectrum 
in KamLAND can be calculated by 

f(E pT ; sin 2 29, Am 2 ) = 

POO 

M / &E V <t(E v ) £ <P 3 (E V )P 3 (E U , sin 2 29, Am 2 ) R(E pi , E' pi ) . (2) 
Jo ■ 

Here R(E pr , E' ) is the energy resolution function and E pT , E' pr are the ob- 
served and the true prompt energies, respectively, and we use a Gaussian 
energy resolution of 7.5%/^£ pr (MeV) Q. The neutrino energy is related to 
the true prompt energy by E v = E' pT + A — m e , where A is the neutron- 
proton mass difference and m e is the positron mass. The cross section o~(E u ) 
for the detection process v e + p — > e + + n is taken from Ref. |l9j |. The neu- 
trino spectrum <p(E u ) from nuclear reactors is well known. I am using the 
phenomenological parameterisation from Ref. |20| and the average fuel com- 
position for the nuclear reactors as given in Ref. The sum over j in Eq. (2) 
runs over 16 nuclear plants, taking into account the different distances from 
the detector Lj and the power output of each reactor (see Table 3 of Ref. 
Finally, Pj(E u , sin 2 29, Am 2 ) is the survival probability for neutrinos emitted 
at the reactor j, depending on the distance Lj, the neutrino energy and the 
two- flavour oscillation parameters sin 2 29 and Am 2 . 

The total number of events predicted for oscillation parameters sin 2 29 and 
Am 2 above the geo-neutrino cut E cut = 2.6 MeV is given by 

POO 

iV pred (sin 2 29, Am 2 ) = / dE pr f{E pT ; sin 2 29, Am 2 ) . (3) 

J -Ecut 

The over-all constant in Eq. (2) is determined by normalising the number 
of events for no oscillations to A^ pre d(sin 2 29 = 0, Am 2 = 0) = 86.8 0. The 
probability distribution of the expected events, i.e. the probability to obtain 
an event with the prompt energy E pr in the interval [E pT , E pr + dE pr ], can be 
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Table 1 



Prompt energies of the 86 anti- neutrino events in KamLAND. 



obtained by normalising the spectrum given in Eq. (2): 



4 „ f(E m ; sin 2 20, Am 2 . 



The prompt energies of the 86 events observed in KamLAND can be extracted 
from Fig. 3 of Ref. 0] and are listed in Tab. 1. Using the 54 events above 
the geo-neutrino cut with E l w > E cnt one obtains the likelihood function 
containing the spectral shape information of the data: 

86 

£ shape (sin 2 29, Am 2 ) = J[ p^; sin 2 29, Am 2 ) . (5) 

i=33 

To take into account also the information implied by the total number of 
observed events I apply the modified likelihood method (see, e.g., Ref. [ttI]): 

£ tot (sin 2 29, Am 2 ) = £ shape (sin 2 29, Am 2 ) x £ ratc (sin 2 29, Am 2 ) (6) 
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with 2 



1 / iV pred (sin 2 29,Am 2 )-N ohs " 

2 \ 0"rate 



£ rate (sin 2 29, Am 2 ) = — exp 

V 2vT(T ra te 

Here iVobs = 54 is the observed number of events, and 

a r 2 ate (sin 2 29, Am 2 ) = iV pred (sin 2 29, Am 2 ) + a 2 yst iV 2 red (sin 2 29, Am 2 ) , (8) 

with the systematical error <T syst = 6.42% 1]. Note that we derive <7 rate from 
the predicted number of events, which introduces the parameter dependence 

Of CTrato- 

By maximising the likelihood function Eq. (6) the best fit parameters Am 2 = 
7.05 x 10~ 5 eV 2 and sin 2 29 = 0.98 are obtained, in very good agreement with 
the values obtained by the KamLAND collaboration: Am 2 = 6.9 x 10 _5 eV 2 
and sin 2 29 = 1 [1]. To calculate allowed regions for the parameters by means 



of the x -cut method one defines 



17, 1 



AX 2 (sin 2 2#, Am 2 ) = 21n£ totimax - 21n£ tot (sin 2 29, Am 2 ) , (9) 

where Aot,max is the maximum of the likelihood function with respect to sin 2 29 
and Am 2 . The 95% confidence regions obtained from Eq. (9) according to 
Eq. (1) are shown in Fig. 1. One finds that they are in excellent agreement 
with the ones published by the KamLAND collaboration. 



2.2 Least-squares analysis of KamLAND data 



The least-squares analyses performed by many authors 0, 0, H, H, 0, HL IE 
loL 11 , 12l | are based on the KamLAND data binned into 13 energy intervals, 
1, . . . , 13, as given in Fig. 5 of Ref. 0]. 3 Since the number of events 



N. 



obs; 



in the individual bins is rather small (in some bins even zero) the use of a 



2 In general a Poisson distribution has to be used for £ ra t e - However, for a mean of 
order iV bs = 54 the Poisson distribution is very well approximated by the Gaussian 
distribution. 

3 In Ref. 46] a different likelihood analysis of KamLAND data has been presented, 
including a goodness of fit evaluation by Monte Carlo methods. Note however, 
that in Ref. [16( the likelihood function is also calculated from the energy binned 
data similar to the least-squares method, in contrast to the event-by-event based 
likelihood discussed in the present work. 
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Figure 1. Comparison of the 95% CL regions obtained from the likelihood and the 
least-squares methods (lines) with the ones published by the KamLAND collabora- 
tion (regions labeled as "Rate+Shape allowed"). The best fit points are marked 
by a star, triangle, dot for the likelihood, least-squares, KamLAND analyses, re- 
spectively. 

least-squares statistic based on the Poisson distribution is appropriate [l8| : 

X 2 (sin 2 29, Am 2 ) = 2 £ |«iV; red - JV*. + iV* bs In -^-| + , 

(10) 

where the term containing the logarithm is absent in bins with no events. The 
predicted number of events N^ red (sm 2 29 , Am 2 ) in bin i is obtained by inte- 
grating the spectrum Eq. (2) over the prompt energy interval corresponding to 
that bin. Eq. (10) has to be minimised with respect to a in order to take into 
account the overall uncertainty of the theoretical predictions <T syst = 6.42% 
Although the "least-squares" character of the statistic X 2 in Eq. (10) is not 
explicitly visible due to the use of the Poisson distribution it is denoted here 
by this term to stress the analogy to the commonly used "x 2 -f u nction" • A 
comparison of KamLAND analvses using Gaussian and Poisson least-squares 
functions can be found in Ref. 

The best fit parameters obtained by minimising Eq. (10) are Am 2 = 7.24 x 
10~ 5 eV 2 and sin 2 20 = 0.90. Assuming that 

AX 2 (sin 2 29, Am 2 ) = X 2 (sin 2 29, Am 2 ) - X^ in (11) 

follows a ^-distribution with 2 degrees of freedom approximate confidence 
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regions are obtained by considering contours of constant AX 2 according to 
the x 2_cu t method in Eq. (1). From Fig. 1 we find that the 95% confidence 
regions obtained by this method are significantly larger than the ones from the 
likelihood analysis and the regions published by the KamLAND collaboration. 
This fact was already noted in Ref. I conclude that the loss of information 
implied by the binning of the data is not negligible, and the likelihood analysis 
provides a more powerful method to extract information from the current 
KamLAND data sample. Note, however, that in future the differences between 
the two methods are expected to decrease, since if more data is available a 
smaller bin size can be chosen, and in the limit of zero bin width the least- 
squares method converges to the un-binned likelihood method. 



3 Confidence regions with correct coverage 

Since the number of events in the currently available KamLAND data sam- 
ple is rather small the question arises, whether the standard procedures to 
calculate confidence regions as described in Sec. 2 are reliable. Especially the 
assumption concerning the distribution of AX 2 might be not justified. More- 
over, the parameters of interest, Am 2 and sin 2 29, enter the problem in a highly 
non-linear way, which leads to multiple local maxima of the likelihood func- 
tion (or local minima of the X 2 -statistic). In such a case the actual confidence 
level of the parameter regions obtained from Eq. (1) may differ significantly 
from the canonical value (3. To check the robustness of the results I have calcu- 
lated frequentist confidence regions, where the correct coverage is guaranteed 
by construction. To this aim I follow the prescription given by Feldman and 
Cousins in Ref. Q. 

For both approaches discussed above - likelihood as well as least-squares meth- 
ods - many synthetic data sets are simulated for fixed oscillation parameters. 
In the case of the likelihood method first the number of events in the synthetic 
data sample iV S i m (sm 2 29, Am 2 ) is generated from a Gaussian distribution with 
mean N pred (sm 2 29 , Am 2 ) and standard deviation cx rate given in Eq. (8). Then 
the prompt energies of the N sim events is thrown according to the distribu- 
tion Eq. (4). To test the least-squares method a value a S i m for the parameter 
a describing the normalisation uncertainty in Eq. (10) is generated from a 
Gaussian distribution with mean 1 and standard deviation cr syst . Then the 
number of events in each bin i is simulated from a Poisson distribution with 
the mean a sim iVp red (sin 2 29, Am 2 ). 

Each "data set" generated this way is analysed as described in Sec. 2 in order 
to calculate AX 2 (sin 2 29, Am 2 ). For each point on a sufficiently dense grid 
in the (sin 2 29, Am 2 ) plane this has been done 10 4 times for the likelihood 
method and 10 5 times for the least-squares method to map out the actual 
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Figure 2. Comparison of confidence regions with correct coverage (lines) with the 
regions obtained from the y^-cut approximation (shaded regions) for the likelihood 
method. 

distribution of AX 2 in that point: p S i m (AX 2 ; sin 2 29, Am 2 ). Then, in analogy 
to Eq. (1), the point (sin 2 29, Am 2 ) is included in the confidence region at (3 
CL if AX| ata (sin 2 29, Am 2 ) obtained from the real data is smaller than the 
one of 100/?% of the simulated data sets in that point in the parameter space: 



The results of these analyses are shown in Figs. 2 and 3 for the likelihood and 
the least-squares methods, respectively. The regions with correct coverage are 
compared to the ones obtained by the ^-cut approximation. In both cases rea- 
sonable agreement of the exact and approximate confidence regions is found, 
although quantitative differences are visible. For the likelihood method the re- 
gions at 99.73% and 99% CL are in excellent agreement, whereas for lower con- 
fidence levels the ^-cvX approximation gives regions somewhat smaller than 
the exact ones. Especially the region 10~ 5 eV 2 < Am 2 < 4 x 10~ 5 eV 2 does 
not appear at 90% CL for the x 2_cu t approximation. In the case of the least- 
squares method the 95% CL regions are in excellent agreement. For higher 
confidence levels the ^-cut approximation gives regions somewhat larger than 
the exact ones, whereas for lower confidence levels the allowed regions are a bit 
underestimated. E.g., the region Am 2 > 2x 10~ 4 eV 2 does not appear at 68.3% 
CL for the ^-cvX approximation. In general the ^-cvX approximation works 
quite well in the vicinity of the best fit point around Am 2 « 7 x 10 -5 eV 2 . 




(12) 



A better understanding of these results can be obtained by considering how 
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Figure 3. Comparison of confidence regions with correct coverage (lines) with the re- 
gions obtained from the x 2_cu t approximation (shaded regions) for the least-squares 
method. 
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Figure 4. Contours of constant as defined in Eq. (12) for the least-squares method 
for (5 = 0.95 (left panel) and (3 = 0.683 (right panel). The thick lines correspond 
to the values which would result from a ^-distribution with 2 degrees of freedom. 
The shaded regions are the approximate allowed regions from the x^-cnt method. 

Xj defined in Eq. (12) varies as a function of the oscillation parameters. Note 
that this quantity does not depend on the actually observed data; it char- 
acterises the properties of the statistical method applied to the specific ex- 
perimental setup. For definiteness I consider the least-squares method, 4 for 
which contours of constant Xj are shown in Fig. 4 for 68.3% and 95% CL. 
In the left panel of that figure one can see that the contour for Xq 95 = 5.99, 



4 Similar behaviour is also found for the likelihood method. 
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which corresponds to the ^-distribution for 95% CL, happens to be rather 
close to the 95% CL region from the ^-cut approximation. This explains the 
good agreement observed in the middle panel of Fig. 3. Furthermore, one finds 
from Fig. 4 that A^ decreases for small values of sin 2 29 and Am 2 . The rea- 
son for this behaviour can be understood as follows: If sin 2 29 becomes < 0.2 
and/or Am 2 < 10 -5 eV 2 the effect of oscillations gets very small, and the sig- 
nal in KamLAND corresponds roughly to the no-oscillation case. Analysing 
data generated from parameters in that region leads to best fit points also in 
the no-oscillation region, with a rather similar X 2 . Hence, the distribution of 
A A 2 is more peaked at low values, which implies relatively small values of Aj. 
This explains the tendency of more constraining exact confidence regions for 
small values of sin 2 29 and Am 2 . The physical reason for this behaviour is that 
even with the present KamLAND data sample no-disappearance can be very 
well distinguished from oscillations with sin 2 29 > 0.2 and Am 2 > 10~ 5 eV 2 . 
Moreover, in the region of small sin 2 29 or Am 2 the full 2-parameter depen- 
dence of the survival probability is lost, and the effective number of degrees of 
freedom is reduced, leading to smaller values of Al. The reason why for the 
likelihood method the agreement is better for higher CL than for the least- 
squares method can be partially attributed to the fact that for the latter the 
approximate regions extend to smaller values of sin 2 29, which implies smaller 
values of Aj and larger disagreement with the ^-approximation. 

In contrast, for sin 2 29 > 0.4 and Am 2 > 10~ 5 eV 2 the simulation yields rela- 
tively higher values of Aj. In that region oscillations are important. However, 
because of the rather small data sample the signature of given parameters can 
not be identified with sufficiently high significance. This implies that statisti- 
cal fluctuations can lead easily to best fit points in a different local minimum. 
In other words, when data are generated by given parameters in that region, 
fluctuations can mimic a signal which is better fitted by quite different param- 
eters. Hence, the best fit points are stronger affected by fluctuations, resulting 
into a broader distribution of A A 2 and larger values of A 2 . From Fig. 4 one 
observes that the approximate confidence regions at 68.3% CL are well inside 
the high Aj regime. This explains the weaker constraints from exact allowed 
regions for low confidence levels in Fig. 3. Note, however, that in the region 
around Am 2 ~ 7 x 10~ 5 eV 2 , where KamLAND is most sensitive to oscilla- 
tions, the A 2 values decrease again. This shows that in that region KamLAND 
can well identify the parameters, leading already to statistical properties close 
to the expected ^-distribution. Therefore, around the best fit point exact and 
approximate confidence regions are in good agreement. 

The main conclusion to be drawn from Figs. 2 and 3 is that already using 
the current 54 events from KamLAND the approximate ^-cnt method gives 
a rather reliable determination of allowed regions. The small quantitative dif- 
ferences to exact confidence regions can be understood by considering the 
impact of statistical fluctuations on the fit. One may expect that once more 
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data will have been collected the differences will further decrease, because 
statistical fluctuations will be less important. Furthermore, if the true oscil- 
lation parameters happen to be close to the present best fit region around 
Am 2 ~ 7 x 10~ 5 eV 2 and large sin 2 29 a rather clear oscillation signal can be 
observed in KamLAND. In that case the impact of other local minima will 
become small and one may expect that AX 2 will follow a distribution rather 
close to the ^-distribution. 



4 Have we already observed oscillations in KamLAND? 

Because of the limited statistics of the current KamLAND data sample the 
data is consistent with an energy independent suppression of the reactor neu- 
trino flux This is evident, since allowed regions appear for large Am 2 , 
corresponding to energy averaged oscillations. However, even this first data 
set indicates some spectral distortion which is consistent with neutrino oscil- 
lations. In this section the statistical significance for an oscillatory signal is 
quantified by using a so-called decoherence parameter. The survival proba- 
bility for the electron anti-neutrinos is modified (in a rather arbitrary way) 
by multiplying the quantum mechanical interference term which leads to the 
oscillatory behaviour by a factor (1 — £): 

P = 1 - - sin 2 29 
2 

Restricting ( to the interval [0, 1] one can describe in a model independent way 
a loss of quantum mechanical coherence due to some unspecified mechanism. 
C = corresponds to standard quantum mechanics, whereas £ = 1 describes 
complete decoherence, i.e., an energy and baseline independent suppression of 
the flux. 5 

Now the data is analysed as a function of the three parameters Am 2 , sin 2 29 
and £. The A A 2 marginalised with respect to Am 2 and sin 2 29 is shown in 
Fig. 5 for the likelihood and the least-squares method. A clear indication in 
favour of oscillations is observed. Complete decoherence is disfavoured with 
AX 2 = 4.4 using the likelihood method and AX 2 = 3.2 by the least-squares 
method. This confirms that the likelihood method is a more powerful tool to 
extract spectral information from the data, in agreement with the results of 
Sec. 2. Assuming that AX 2 (() is distributed as a x 2 with 1 degree of freedom 

5 Decoherence might have its origin e.g. in quantum gravity See also Ref. [23| 
and references therein, for an application to atmospheric neutrino oscillations. A 
method similar to Eq. (13) has been used in Ref. |2J] to investigate the evidence for 
quantum mechanical interference in the B°B and K°K systems. 



1 - (1 - C) cos 



Am 2 L 



2K 



(13) 
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Figure 5. AX 2 as a function of the decoherence parameter £ for the likelihood and 
the least-squares method. 

I conclude that the current KamLAND data provides a ~ 2a indication in 
favour of neutrino oscillations, 6 implying quantum mechanical interference 
over distances of the order of 200 km. 

Depending on the true values of the oscillation parameters one may expect 
that the statistical significance for oscillations in KamLAND will strongly 
improve by future data. A simple rescaling of the current data by a factor 5 
leads to an exclusion of complete decoherence at 3.7a if the true value of Am 2 
turns out to be 7 x 10 -5 eV 2 . On the other hand if Am 2 = 1.5 x 10" 4 eV 2 
decoherence can be excluded only at 2.6a, since for large Am 2 the baselines 
in KamLAND are too long to be sensitive to the oscillations. 



5 Conclusions 



In this letter two different analysis methods for the first data from the Kam- 
LAND reactor neutrino experiment have been compared. I found that an 
event-by-event based likelihood method provides a more powerful tool to ex- 
tract information on two-neutrino oscillation parameters than a least-squares 
method based on energy binned data. The likelihood method takes into ac- 
count the precise energy information contained in each single event and avoids 
the information loss due to binning. Furthermore, exact frequentist confidence 

6 Note that here a relative comparison of the fits for oscillations and decoherence 
is performed; no statement about the absolute qualitiy of the fit is made. Hence, 
these results are in agreement with Ref. 0], where the observed spectrum is found 
to be consistent with an oscillation signal at 93% CL but with flat suppression at 
53% CL. 
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regions in the parameter space have been calculated by means of Monte Carlo 
simulation according to the Feldman- Cousins method [14]. This method prop- 
erly accounts for the non-linearity of the oscillation parameters sin 2 29 and 
Am 2 , and statistical fluctuations in the data, which can be quite large due to 
the rather small number of events in the current data sample. I have found a 
reasonable agreement of the exact confidence regions with the ones obtained 
from the x 2 -cut approximation, especially in the vicinity of the best point at 
Am 2 fa 7 x 10~ 5 eV 2 . However, depending on the analysis method (likelihood 
or least-squares) quantitative differences are visible, especially for lower confi- 
dence levels and far from the best fit point. Finally, although the current data 
is consistent with an energy independent flux suppression, a ~ 2a indication 
in favour of oscillations can be stated using the likelihood method, which is 
especially sensitive to the spectral shape information. Put in other words, this 
implies a ~ 2a evidence for quantum mechanical interference over distances 
of the order of 200 km. 

In summary, the results obtained in this work confirm that even for the limited 
statistics of the current KamLAND data sample the x 2_cu t approximation to 
calculate confidence regions for the oscillation parameters gives rather reli- 
able results. One expects that in future the agreement between approximate 
and exact confidence regions will improve due to increase in statistics. More- 
over, the differences between likelihood and least-squares methods will become 
smaller. 
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